/* @(#)trand.c 1.6 1/13/95 */
/* trand.c 23-5-96  adapted from ``Numerical Recipes in C'' by Press,
   Flannery, Teukolsky and Vetterling.  Tdrand() returns a uniform random
   deviate between 0.0 and 1.0.  Tlrand() returns a uniform random deviate
   between 0 and 2^31 */
/* MemProf labels: #18XX, in use: none */
/* FEError labels: #18XX, in use: 1800 */

#include "tierra.h"
#include "globals.h"

#define M1 259200L
#define IA1 7141L
#define IC1 54773L
#define RM1 (1.0/M1)
#define M2 134456L
#define IA2 8121L
#define IC2 28411L
#define RM2 (1.0/M2)
#define M3 243000L
#define IA3 4561L
#define IC3 51349L

/*
 * n a m e
 *
 * description
 *
 * Returns:
 */
void tsrand(seed)
    I32s seed;
{   I32s j;

    RandIx1 = (IC1 + seed) % M1;  /* Seed the first routine */
    RandIx1 = (IA1 * RandIx1 + IC1) % M1;
    RandIx2 = RandIx1 % M2;  /* and use it to seed the second */
    RandIx1 = (IA1 * RandIx1 + IC1) % M1;
    RandIx3 = RandIx1 % M3;  /* and third routines. */
#ifdef __MWERKS__
    if (RandIx1 < 0)
        RandIx1 *= -1;
    if (RandIx2 < 0)
        RandIx2 *= -1;
    if (RandIx3 < 0)
        RandIx3 *= -1;
#endif /* __MWERKS__ */
    for (j = 1; j <= 97; j++)  /* Fill the table with sequential uniform */
    {   RandIx1 = (IA1 * RandIx1 + IC1) % M1;  /* deviates generated by the */
        RandIx2 = (IA2 * RandIx2 + IC2) % M2;  /* first two routines. */
        TrandArray[j] = (RandIx1 + RandIx2 * RM2) * RM1;
    }           /* Low- & high-order pieces combined here */
}

/*
 * n a m e
 *
 * description
 *
 * Returns:
 */
double tdrand()
{   double temp;
    I32s j;

    RandIx1 = (IA1 * RandIx1 + IC1) % M1;
    /* Except when initializing, this is where we */
    RandIx2 = (IA2 * RandIx2 + IC2) % M2;
    /* start.  Generate the next number for each */
    RandIx3 = (IA3 * RandIx3 + IC3) % M3;  /* sequence. */
    j=(I32s)(1+((97*RandIx3)/M3));
    /* Use the third sequence to get an integer between */
    if (j > 97 || j < 1)  /* 1 and 97 */
        FEError(-1800, EXIT, WRITE, "Tierra tdrand() This cannot happen.");
    temp = TrandArray[j];  /* Return that table entry, */
    TrandArray[j] = (RandIx1 + RandIx2 * RM2) * RM1;  /* and refill it. */
    return temp;
}

#ifdef FUTURE

/*
 * n a m e
 *
 * description
 *
 * Returns:
 */
float gasdev()
{   static I32s iset = 0;
    static float gset;
    float fac, r, v1, v2;

    if (iset == 0)
    {   do
        {   v1 = 2.0 * tdrand() - 1.0;
            v2 = 2.0 * tdrand() - 1.0;
            r = v1 * v1 + v2 * v2;
        }
        while (r >= 1.0);
        fac = sqrt(-2.0 * log(r) / r);
        gset = v1 * fac;
        iset = 1;
        return v2 * fac;
    }
    else
    {   iset = 0;
        return gset;
    }
}
#endif /* FUTURE */

#undef M1
#undef IA1
#undef IC1
#undef RM1
#undef M2
#undef IA2
#undef IC2
#undef RM2
#undef M3
#undef IA3
#undef IC3
